% -------------------------------------------------------------------------
% Function: [class,type]=slowdbscan(x,k,Eps)
% -------------------------------------------------------------------------
% Aim: 
% Clustering the data with Density-Based Scan Algorithm with Noise (DBSCAN)
% -------------------------------------------------------------------------
% Input: 
% x - data set (m,n); m-objects, n-variables
% k - number of objects in a neighborhood of an object 
% (minimal number of objects considered as a cluster)
% Eps - neighborhood radius, should be a vector(n*1), which indicates the 
% radius at dim x in [1,n].
% -------------------------------------------------------------------------
% Output: 
% class - vector specifying assignment of the i-th object to certain 
% cluster (m,1), 0 indicates cluster, >=1 indicates a class.
% -------------------------------------------------------------------------
% Example of use:
% x=[randn(30,2)*.4;randn(40,2)*.5+ones(40,1)*[4 4]];
% [class,type]=slowdbscan(x,5,[])
% clusteringfigs('slowdbscan',x,[1 2],class,type)
% -------------------------------------------------------------------------
% References:
% [1] M. Ester, H. Kriegel, J. Sander, X. Xu, A density-based algorithm for 
% discovering clusters in large spatial databases with noise, proc. 
% 2nd Int. Conf. on Knowledge Discovery and Data Mining, Portland, OR, 1996, 
% p. 226, available from: 
% www.dbs.informatik.uni-muenchen.de/cgi-bin/papers?query=--CO
% [2] M. Daszykowski, B. Walczak, D. L. Massart, Looking for 
% Natural Patterns in Data. Part 1: Density Based Approach, 
% Chemom. Intell. Lab. Syst. 56 (2001) 83-92 
% -------------------------------------------------------------------------
% Written by Michal Daszykowski
% Department of Chemometrics, Institute of Chemistry, 
% The University of Silesia
% December 2004
% http://www.chemometria.us.edu.pl

% cloned from https://github.com/sinjax/dbscan
% modified by xuyuntao at 20230115

function [ptClass] = slowdbscan_dist(pts,minPtNums,maxDis)
% maxDis为n*1向量，表示各维距离

    [ptNums,dims] = size(pts);   % 点数和维数

    if not(length(maxDis) == dims)
        error('The dimension of maxDis should be equal to or larger than size(pts,2).');
    end

    maxDis = reshape(maxDis,1,[]);   % 将maxDis转为1行，方便计算

    ptClsTmp = 1;   % 当前点的聚类索引
    touched = false(ptNums,1);   % 该点是否已访问
    ptClass = zeros(ptNums,1);   % 各点聚类索引
    ptIdxCheck = false(ptNums,1);  % 该聚类待检查点索引，对应pts行索引
    
    for ptIdx = 1:ptNums
        ptIdxCheck(:,1) = false;
        if not(touched(ptIdx,1))  % 只有未访问的点才进行计算
            ob = pts(ptIdx,:);
            dis = abs(bsxfun(@minus,ob(1,:),pts(:,:)));
            checkPtsIdx = all(bsxfun(@le,dis,maxDis),2);
%             checkPtsIdx = (dis <= maxDis);
            ptIdxCheck = or(ptIdxCheck,checkPtsIdx);   % 动态维护聚类该类的待检查点
            
            if nnz(ptIdxCheck) >= minPtNums   % 大于最小类点数的聚类
                ptClass(ptIdxCheck,1) = ptClsTmp;  % 该点为一个聚类，标记其聚类索引
                
                while any(ptIdxCheck)
                    ptIdxCheckFirst = find(ptIdxCheck,1,'first');

                    obSub = pts(ptIdxCheckFirst,:);
                    touched(ptIdxCheckFirst,1) = true;

                    ptIdxCheck(ptIdxCheckFirst,1) = false;

                    D = abs(bsxfun(@minus,obSub(1,:),pts(:,:)));
                    ptIdxCheckSub = all(bsxfun(@le,D,maxDis),2);
%                     ptIdxCheckSub = (D <= maxDis);
                    
                    if any(ptIdxCheckSub)
                        ptClass(ptIdxCheckSub,1) = ptClsTmp;
                        
                        ptIdxCheckSubValid = and(not(touched),ptIdxCheckSub);

                        touched(ptIdxCheckSubValid,1) = true;
                        ptClass(ptIdxCheckSubValid,1) = ptClsTmp;

                        ptIdxCheck(ptIdxCheckSubValid,1) = true;
                    end
                end
                ptClsTmp = ptClsTmp + 1; 
            end
        end
    end

end


%...........................................
function [Eps]=epsilon(x,k)

% Function: [Eps]=epsilon(x,k)
%
% Aim: 
% Analytical way of estimating neighborhood radius for DBSCAN
%
% Input: 
% x - data matrix (m,n); m-objects, n-variables
% k - number of objects in a neighborhood of an object
% (minimal number of objects considered as a cluster)

    [m,n]=size(x);
    
    Eps=((prod(max(x)-min(x))*k*gamma(.5*n+1))/(m*sqrt(pi.^n))).^(1/n);

end


%............................................
function [D]=dist(i,x)

% function: [D]=dist(i,x)
%
% Aim: 
% Calculates the Euclidean distances between the i-th object and all objects in x	 
%								    
% Input: 
% i - an object (1,n)
% x - data matrix (m,n); m-objects, n-variables	    
%                                                                 
% Output: 
% D - Euclidean distance (m,1)

    [m,n]=size(x);
    D=sqrt( sum( (((ones(m,1)*i)-x).^2) , 2) );
    
    if n==1
       D=abs((ones(m,1)*i-x))';
    end

end